input=imread('Cap_50_4.2501.bmp','bmp');
%input=input'/255;

lz=75;
lx=200;
ly=50;
res=92;
startp=3750;

ino=1-input(startp:startp+lx-1,1:ly,2)/255;

figure
imagesc (ino); figure(gcf);

figure

a=zeros(lx/2,ly/2);

for i=1:2:lx-1
    for j=1:2:ly-1
        if (ino(i,j)+ino(i+1,j)+ino(i,j+1)+ino(i+1,j+1)>2)
            a((i+1)/2,(j+1)/2)=1;
        else
            a((i+1)/2,(j+1)/2)=0;
        end
    end
end

imagesc (a); figure(gcf);

aa=zeros(lx/2+res,ly/2+2);

for i=1:lx/2
    for j=1:ly/2+2
        if ((j==1) || (j==ly/2+2))
            aa(i,j)=1;
        else
            aa(i,j)=a(i,j-1);
        end
        
    end
end





 fid = fopen('rough_capillary.vtk','wt');
% fprintf(fid,'# vtk DataFile Version 2.0\n');
% fprintf(fid,'J.Yang Lattice Boltzmann Simulation 3D Capillary Trapping\n');
% fprintf(fid,'ASCII\n');
% fprintf(fid,'DATASET STRUCTURED_POINTS\n');
% fprintf(fid,'DIMENSIONS         %i         %i         %i\n',lx/2+res,ly/2+2,lz);
% fprintf(fid,'ORIGIN 0 0 0\n');
% fprintf(fid,'SPACING 1 1 1\n');
% fprintf(fid,'POINT_DATA     %i\n',(lx/2+res)*(ly/2+2)*lz);
% fprintf(fid,'SCALARS sample_scalars float\n');
% fprintf(fid,'LOOKUP_TABLE default\n');     

 for k=1:lz
     for j=1:ly/2+2
         for i=1:lx/2+res
             
         if ((k==1) ||(k==lz)) && (i<=lx/2)
             fprintf(fid,'%1d\n',1);
         else
             fprintf(fid,'%1d\n',aa(i,j));
         end
         
         end
     end
 end
 
 lx/2+res
 ly/2+2
 lz
